// 2D SPDE model of thin film deposition model. See Ni Dong's dissertation for more details.
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip>
using namespace std;
#include "model_SPDE2D.h"

const int ORDER = 2;
// The following constant is used in update()
const double k_w = 3.3829e-12;
const double a_w = 0.6042;
const double E_w = 2.7e-3; // Unit = EV
const double k_c = 1.0274e-13;
const double a_c = 0.1669;
const double E_c = 1.9e-3; // Unit = EV
const double a_v = 15.55493;
const double k_v = 20.64504;
const double a_t = 0.02332;
const double k_t = 0.0261;
const double k_B = 8.617343e-5; // Boltzmann's constant, Unit = ev/K
// ---------------------------------------------------
model_SPDE2D::model_SPDE2D(int iMode,int iLatticeSize,double idt)
{
	int i_dim[2] = {1,1};
	int o_dim[1] = {iMode*iMode};
	block::initialize(2,i_dim,1,o_dim,idt);

	int i,m,n;
	mode = iMode;
	LatticeSize = iLatticeSize;
	c = new double[ORDER];
	sigma2 = 0.0;
	lambda_Re = new double*[mode];
	lambda_Im = new double*[mode];
	Km2 = new double*[mode];
	for(i=0;i<mode;i++){
		lambda_Re[i] = new double[mode];
		lambda_Im[i] = new double[mode];
		Km2[i]       = new double[mode];
	}
	//z_mean = new double[mode*mode];
	z_mean = new double[mode*mode];
	z_cov  = outputs_next[0];

	double dx = pi/LatticeSize;
	double temp = dx*dx*pi*pi;
	for(m=0;m<mode;m++){
		for(n=0;n<mode;n++){
			lambda_Re[m][n] = 0.0;
			lambda_Im[m][n] = 0.0;
			z_mean[index2D(m,n)] = 0.0;
			z_cov[index2D(m,n)]  = 0.0;
			Km2[m][n] = (2-2*cos(2*m*dx))/temp;
		}
	}
}
// ---------------------------------------------------
model_SPDE2D::~model_SPDE2D()
{
	int i;
	
	delete [] z_mean;
	//delete [] z_cov;

	for(i=0;i<mode;i++){
		delete [] lambda_Re[i];
		delete [] lambda_Im[i];
	}
	delete [] lambda_Re;
	delete [] lambda_Im;
	delete [] c;
}
// ---------------------------------------------------
void model_SPDE2D::update(double sysTime, double W, double T,const double *z_mean0)
{	
	int m,n;

	if(z_mean0 != NULL){
		for(m=0;m<mode;m++){
			for(n=0;n<mode;n++){
				z_mean[index2D(m,n)] = z_mean0[index2D(m,n)];
			}
		}
	}	

	while(time<=sysTime)
	{
		block::output();
		
		c[0] = W*(1-k_w/(pow(W,a_w)*exp(-k_B*T/E_w))); //eq(5.28)
		c[1] = k_c/(LatticeSize*LatticeSize*pow(W,a_c)*exp(-k_B*T/E_c));   //eq(5.29)
		sigma2 = pi*pi/(LatticeSize*LatticeSize)*(1+exp((a_t+k_t*W)*T)/exp(a_v+k_v*W)); //eq(5.30)

		for(m=0;m<mode;m++){
			for(n=0;n<mode;n++){
				lambda_Re[m][n] = -4*(m*m+n*n)*c[1]; //eq(5.10)
				if(m!=0||n!=0)
				{
					z_mean[index2D(m,n)] = exp(lambda_Re[m][n]*dt)*z_mean[index2D(m,n)];//eq(5.16)
					z_cov[index2D(m,n)] = exp(2*lambda_Re[m][n]*dt)*z_cov[index2D(m,n)]+sigma2*(exp(2*lambda_Re[m][n]*dt)-1)/(2*lambda_Re[m][n]);//Modified eq(5.16) and eq(5.17)
					//z_cov[index2D(m,n)] = sigma2*(exp(2*lambda_Re[m][n])-1)/(2*lambda_Re[m][n])-z_mean_old*z_mean_old+z_mean[0][0]*z_mean[0][0];//eq(5.16) and eq(5.17)
					// z_mean[index2D(m,n)] = 0;
					// z_cov[index2D(m,n)] = sigma2*(exp(2*lambda_Re[m][n]*(time+dt))-1)/(2*lambda_Re[m][n]);
				}
				else
				{
					z_mean[index2D(0,0)] = z_mean[index2D(0,0)]+dt*pi*c[0];
					// z_cov[index2D(0,0)] = z_cov[index2D(0,0)]+sigma2*dt;
					z_cov[index2D(0,0)] = z_mean[index2D(0,0)]*z_mean[index2D(0,0)];
					//z_mean[0][0] = (time+dt)*pi*c[0];
					//z_cov[0][0] = sigma2*(time+dt);
				}
			}
		}	
		
		time += dt;
	}
}
// ---------------------------------------------------
void model_SPDE2D::update(double sysTime)
{
	double W = inputs[0][0];
	double T = inputs[1][0];
	update(sysTime,W,T);
}
// ---------------------------------------------------
double model_SPDE2D::getR2()
{
	int m,n;
	
	double r2 = 0.0;
	for(m=0;m<mode;m++){
		for(n=0;n<mode;n++){
			if(m!=0||n!=0){
				r2+=z_cov[index2D(m,n)];
			}
		}
	}
	r2 = r2/(pi*pi); // eq(5.34)
	return r2;
}
// ---------------------------------------------------
double model_SPDE2D::getM2()
{
	int m,n;
	double m2 = 0.0;
	for(m=0;m<mode;m++){
		for(n=0;n<mode;n++){
			if(m!=0||n!=0){
				m2+=Km2[m][n]*z_cov[index2D(m,n)];
			}
		}
	}
	return m2;
}
// ---------------------------------------------------
double model_SPDE2D::getH()
{
	return z_mean[index2D(0,0)]/pi;
}
// ---------------------------------------------------
void model_SPDE2D::getZmean(const string filename)
{
	ofstream fout(filename.c_str());
	fout << "mode = " << mode << "\n";
	int m,n;
	for(m=0;m<mode;m++){
		for(n=0;n<mode;n++){
			fout << setw(10) << z_mean[index2D(m,n)];
		}
		if(m!=mode-1)
			fout << "\n";
	}
	fout.close();
}
// ---------------------------------------------------
void model_SPDE2D::getZmean(const double* z)
{
	z = z_mean;
}
// ---------------------------------------------------
